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We investigate the structure of the wind in the neutron star X-ray binary system Vela X-1 by analyzing its flaring behavior Vela X-1 
shows constant flaring, with some flares reaching fluxes of more than 3.0 Crab between 20-60 keV for several 100 seconds, while the 
average flux is around 250mCrab. We analyzed all archival INTEGRAL data, calculating the brightness distribution in the 20-60 keV 
band, which, as we show, closely follows a log-normal distribution. Orbital resolved analysis shows that the structure is strongly 
variable, explainable by shocks and a fluctuating accretion wake. Analysis of RXTE ASM data suggests a strong orbital change of 
A'h- Accreted clump masses derived from the INTEGRAL data are on the order of 5 x 10"-10"' g. We show that the lightcurve can 
be described with a model of multiplicative random numbers. In the course of the simulation we calculate the power spectral density 
of the system in the 20-100 keV energy band and show that it follows a red-noise power law. We suggest that a mixture of a clumpy 
wind, shocks, and turbulence can explain the measured mass distribution. As the recently discovered class of supergiant fast X-ray 
transients (SFXT) seems to show the same parameters for the wind, the link between persistent HMXB like Vela X- 1 and SFXT is 
further strengthened. 

Key words. Accretion - X-rays: binaries - X-rays: individual (Vela X-1) - Methods: statistical 
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1. Introduction 

The high-mass X-ray binary (HMXB) Vela X-1 was discovered 
in 1967 ( jChodil et al.| 11967) 1 and is today one of the best studied 
objects of its class hosting a neutron star. It is a bright, eclips- 
ing source, showing strong pulsations with a pulse period of 
around 283 sec ( McClintock et al.||1976[l and an orbit al period of 
8.9 days ( ,van Kerkwijk et al.||1995| l. |Quaintrell et al.| ( i2003 ) have 
shown that the separation between the neutron star and its optical 
companion, the BO.SIb supergiant HD 77581, is only 1.7 stellar 
radii. The neutron star is thus deeply embedded in the stellar 
wind of the s upergiant, which h as a mass loss rate of about 



10 ''Mgyr ' ( [Nagase et al.[|1986^ . Material from the wind is 



accreted by the neutron star and is channeled by its strong mag- 
netic field onto the magnetic poles. The accretion rate leads to 
an average X-ray luminosity of ~ 4 x 10""' erg s"', although the 
luminosity is strongly variable on all time-scales, varying up to 
a factor of a t least 20-30 ([Staubert et al.| |2004[ [Kreykenbohm 
et al.||2008| l. |Staubert et ai| ( |2004| ) found short, but very bright 
flares reaching up to almost 7 Crab (1000 counts sec"' [cps] in 
the 20^0 keV band of INTEGRAL ISGRI). 

The physical processes leading to this strong variation are 
unclear, however. Because the luminosity of the neutron star 
is proportional to the mass accretion rate, changes in the den- 
sity or the effective velocity of the stellar wind directly lead 
to variations in the X-ray luminosity ( jBondi & Hoyle| |1944[ 



Davidson & Ostriker 1973 i. Changes of that kind have shown 
up in simulations of the system ( Blondin et al.j [1990 199T} 



Mauche et al. 2008|l. Strong density changes have also been 



proposed to explain varying line strengths in high-resolution X- 
ray spectra (Sako et al. 1999'; "Watanabe et al. 2006 1 or were 



inferred via measurement of emission lines in hot, thin plasma 
at the same time with fluorescence lines in colder, optical thick 
plasma ( |Schulz et al.j |2002[ [Goldstein et ar]|2004| l. 

The goal of this paper is to study the variations in the 
lightcurve and investigate a possible connection between the 
wind structure and the flaring behavior. This paper is structured 
as follows: in Sect.[2[we present our data and analysis method. In 
Sect. [3] we show an orbital phase averaged and orbital phase re- 
solved analysis of the INTEGRAL ISGRI hard X-ray and RXTE 
ASM soft X-ray lightcurve. In order to constrain the physical 
processes leading to the observed variability we present a numer- 
ical method to simulate lightcurves in Sect.[4| For this simulation 
we calculate and model the power spectral density of Vela X-1. 
In Sect. [5] we summarize our results and discuss them in a larger 
context. 

2. Observations and data 

The "International Gamma-Ray Astrophysics Laboratory" 



(INTEGRAL, Winkler et al. 2003 i carries three X-ray instru 



ments: JEM-X ^Lund et al.||2603| l for the soft X-rays, the spec 
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Table 1, Overview of the INTEGRAL data of Vela X- 1 





Revolutions 


MJD 




Exposure 


Block 1 
Block 2 
Block 3 


137-141 
373-383 
433-440 


52970.4 - 
53678.7 - 
53857.1- 


- 52984.9 
-53708.2 
-53879.0 


0.92 Msec 
1.57 Msec 
1.13 Msec 



trometer SPI ( Vedrenne et al. 2003') and the "Imager on Board 
the INTEGRAL Sate llite" (IBIS) for the hard X- and y-rays 
(|Ubertini et al. 2003 1. All these instruments are coded mask in- 



struments ( Renaud et al. 2006 and references therein), allowing 



imaging of the hard X-ray sky. IBIS has a field of view of almost 
30° X 30° and its detector consists of two planes: the bottom 
layer PICsIT, sensitive in the 200 keV-10 MeV range and the top 
layer ISGRI, which is sensitive between 15 keV-1 MeV (Lebrun 
|et al.| |2003| l. Because of the coded mask instruments, observa- 
tions with INTEGRAL are usually carried out in a dithering 
mode, where the pointing of the satellite is shifted slightly ev- 
ery 1800-3600 sec. Each of these pointings makes up one sci- 
ence window (ScW). During slews between ScWs, taking about 
150 sec, no data are gathered. INTEGRAL has provided a large 
fraction of the detailed hard X-ray studies of Vela X-1, produc- 
ing almost 3.6 Msec of good data between 2003 November and 
2006 May, split mainly into three large blocks separated by some 
hundred days. All observations of INTEGRAL that were public 
as of 2009 December 01 were used as the data basis for this 
work, see Table [T] 



3. Statistical analysis 

3.1. Phase averaged analysis 

While data from individual observation blocks have previously 
been analyzed in detail (Kreykenbohm et al 



200m Schanne 



et al. 2007[ l, studying all available data allowed us to carry out a 
deep statistical analysis of the flaring behavior of Vela X-I. Our 
analysis was performed with special regard to the giant flares 
with peak luminosities of more than six times the average lu- 
minosity ( [Staubert et~al . 2004; Kreykenbohm et al.||2008| l, i.e., 
more than 300 cps in ISGRI (-1.3 Crab) between 20-60 keV. It 
was not clear if these giant flares can be explained by the same 
physical process as the other, common, flares in Vela X-1 or if 
they are caused by a different mechanism. 

For this work, only ISGRI data between 20-60 keV were an- 
alyzed, because in this energy range the source is bright and the 
flux is unaffected by photoabsorption, because typical values of 
the equivalent hydrogen column A^h outside the eclipse are be- 
low 3 X 10^^ cm"^ (ICreykenbohm et al.,_1999j . Our analysis ex- 
tends the energy range of prior observations of the stellar wind, 
which used high-resolution spectra in soft X-rays to investigate 
the features of the wind (Watan abe et al.| |2006| and references 
therein). 

We extracted lightcurves from ISGRI using ii_light, a tool 
distributed as part of the Offline Scientific Analysis (OSA) 7.0 
package. The time resolution was chosen to be 283.5 sec to av- 
erage each data-point over one pulse period to eliminate these 
fluctuations. The pulse period is changing in a random-walk be- 
havior on all times scales ( Deeter et al. 1989| l, but no long-term 
trend is evident, so that the chosen time resolution is sufficient 
for the purpose of this investigation. 

As an example for the analyzed data the upper panel of Fig.fT] 
shows the lightcurve of the 20-60 keV band for data from revo- 
lutions 433-440 (Block 3). 
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Fig. 1. Top: Lightcurve in the 20-60 keV band. The vertical 
dashed lines show the start and the end of the eclipse accord- 
ing to the ephemeris of |Kreykenbohm et al. (j2008J, and the 
dash-dotted line shows the respective center of eclipse. Bottom: 
Simulated lightcurve with the same statistical parameters and 
temporal resolution as the observed lightcurve above, but with- 
out eclipses. See text for details. 
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Fig. 2. Histogram of the lightcurve of the 20-60 keV band, 
binned to 256 bins. The solid curve shows the best-fit single 
Gaussian, the dash-dotted one the histogram of the simulated 
lightcurve as described in the text. The dotted histogram to the 
left shows the noise level of the background. 



The lightcurve data were filtered according to orbital phase 



"orb. 



leaving only data of phases between 0.19 



f'orb 



< 0.81, 



for which the system is out of eclipse. These data were then 
binned into 256 count rate bins. The bins are spaced logarith- 
mically between 1 cps and 1000 cps. This binning leads to a his- 
togram of the orbital phase averaged brightness distribution of 
the source, shown in black in Fig. |2] The distribution closely 
resembles a normal distribution in log-space, i.e., a log-normal 
distribution (F iirst et al. , 200 8). In order to quantify the shape of 
the distribution we fitted a Gaussian function in log-space, i.e., 
a log-normal distribution in countrate space. Following the ap- 
proach by Uttley et al. ( 2005) l the uncertainties of the values A^, of 
each histogram data bin were estimated to be y/Wj, i.e., assum- 
ing Poisson statistics. We show that this assumption is justified 
in Sect. |4] The best-fit function, with a ;t'^-value of 190 for 81 
degrees of freedom, is shown as solid line in Fig. |2] where in the 
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Orbital phase 

Fig. 3. Landscape plot of orbital phase resolved histograms of 
the ISGRI 20-60 ke V data. Color-coded is the probability for a 
datapoint to fall into the respective histogram bin. The black line 
shows the median count rate in each phase bin. Uncertainties of 
that value are also plotted, but they are too small to be clearly 
visible in this plot. Note that the same histograms are shown 
twice for clarity. 



fit only bins with more than 20 measurements were taken into 
account (see |Gehrels| |1986| ). A Gaussian function in log-space 
is best characterized by its mean and standard deviation. These 
values represent the median of the log-normal function, which 
will be denoted as (x), and its multiplicative standard deviation, 
&, i.e., 68.3% of all data points fall in the interval [{x} /&, {x) ■ cr] 



(Ahrens 19541, resulting in asymmetric error bars according to 



the skewness of the distribution. Note that consequently & is 
unitless. With these definitions, the best-fit function has a median 
value of {x)f^i - 50.3 cps and a multiplicative standard deviation 
of CTfit =1.8. Compared with that the measured count rates have 
a median flux of (x)nieas - 48.2 cps and a multiplicative standard 
deviation of (Tmeas = 2.0, similar to the fit. 

The largest cause of source-independent variation that affects 
the data are background fluctuations. To obtain an estimate for 
the background noise level, we extracted a lightcurve in the 20- 
60 keV energy band from a region in the sky without any known 
X-ray source {a = 9h 16m 37s, 6 = -44°07'11"). The bright- 
ness distribution of this background is shown with a dotted line 
in Fig. l2] The median value of the background distribution is 
(j^)bkf ~l-8 cps, which is below 99.98% of the source data. Note 
that the brightest data points of the background reach values of 
up to 11 cps, which is in the lower left flank of the source dis- 
tribution. They do not, however, change the overall shape and 
parameters of the distribution, so that a significant background 
influence can be ruled out. 



3.2. Orbital phase resolved analysis 

We have seen in the previous section that deviations from a 
log-normal distribution are visible. These deviations could be 
induced by averaging over all orbital phases, where possibly 



not all phases show the same statistical brightness variations. 
Systematic influences could be due to the complex and turbu- 
lent accretion geometry, consisting of a bow shock, an accretion 
wake, a possible photoionization wake or even a tidal stream 
(see, e.g., [Blondin et ai:] [1990) [T99T] [Mauche et ai;i[2008;). An 



ladie et al. 



accretion wake in Vela X- 1 was already postulated by[ 

( |1975| l, showing up as absorption dips in the lightcurve of the 
Ariel V Sky Survey Experiment. 

To analyze the orbital influences, we extracted orbital phase 
resolved histograms. The data set comprises ~8 orbits overall, 
which allows us to extract histograms with A^orbit = 0.05. We re- 
duced the number of logarithmically spaced bins to 96 to collect 
more than 20 values in most bins to satisfy the requirements for 
X^ statistics. Figureplshows color-coded the probability of find- 
ing a datapoint in a countrate bin for a given histogram, equiv- 
alent to the height of the histogram in Fig. l2] Additionally the 
average countrate in every bin is shown, together with its uncer- 
tainty. As the average countrate can be calculated from a large 
number of data points in the lightcurve, these uncertainties are 
rather small. Significant variations from phase bin to phase bin 
can be seen in the figure. Especially notable is a dip in the av- 
erage countrate around phase ^ 0.3-0.4 and a sudden bright- 



ening to twice the value afterwards. Goldstein et al. (2004) ana- 
lyzed Chandra data taken at cpoih = 0-25 and (poit = 0.5 of one or- 
bit and measured a stronger absorbed spectrum in the 0orb = 0.5 
data. These authors explained the increased absorption by the 
accretion wake being in the light of sight during that episode. 
The dip in Fig.[3]occurs slightly before (p - 0.5. If an accretion 
wake is responsible for this feature, it seems to have shifted to 
earlier orbital phases in our data. 

When analyzing the three data blocks (Table [Til separately 
an orbital modulation is also visible, although with different 
strength in the diff'erent blocks. This indicates a strongly vari- 
able feature rather than an orbitally stable one and suggests that 
the variation seen in Fig.l3]is biased by single orbits probably de- 
viating strongly from the average orbital histogram, maybe due 
to a particularly active or quiet state of the system. 

The influence of different A^h values during the orbit can be 
investigated by analyzing data in softer X-rays. Softer energy 
bands are far more affected by photoabsorption than the ISGRI 
band, because the cross-section of the photoabsorption is declin- 
ing quickly with increasing energy as cr oc E^^. A very good data 
set in the soft X-rays is provided by the All-Sky Monitor (ASM, 
Levine et al. 1996) aboard the "Rossi X-Ray Timing Explorer" 



(RXTE). The ASM performs 90 sec dwells on an irregular ba- 
sis onto different sources, including Vela X-1, and thus provides 
good statistics in the three energy bands between 1.5-lOkeV 
called A (1.5-3 keV), B (3-5 keV), and C (5-lOkeV) band. 

We performed a similar analysis of the ASM data as of the 
ISGRI data, but binning the countrate directly to 256 bins in- 
stead of its logarithm. The linear binning was necessary because 
especially in the eclipse many data points of the lightcurve are 
negative and cannot be presented on a logarithmic scale. It has 
previously been shown that due to the complex accretion geom- 
etry of the system, the absorption coefficient frequently changes 
dramaticall y over the orbit (Hab erl & Whitej |I990[ [KretschmE] 
let al.||2008|). Especially in the late orbital phases after <potb = 0.5, 
the accretion and photoionization wake, consisting of denser ma- 
terial, are expected to be in the line of sight ( |Kaper et al.||I994) . 
To analyze the influence of the absorption on the brightness dis- 
tribution, we divided the data into the same 20 orbital phases as 
the ISGRI data. 

Figure|4]shows the histograms for the A+B band and for the 
C band. We added the narrow A and B band as both showed very 
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similar behavior, but only low significance. By adding A and B 
the countrates are comparable to the C band. It is obvious that 
in both bands the mean luminosity is monotonically decreasing 
over the orbit (see also We n et al.||2006| l. At the same time, the 
distribution narrows, meaning that strong deviations from the 
mean are more common in the early orbital phases than in the 
later ones. It is remarkable that no sudden dip structures like in 
the ISGRI data in Fig. l3] are visible. This could be because the 
ASM data stretch over more than 13 years, averaging out any in- 
dividual orbit variations. As the ISGRI data comprise only eight 
orbits, they are strongly influenced by single outliers and do not 
provide a sufficient time basis for determining an average orbital 
profile. Nonetheless, no declining trend in the average brightness 
is evident in ISGRI. 

The missing trend strongly points to the explanation that the 
trend in the ASM is due to increasing absorption by a cold ab- 
sorber This interpretation is supported by the analysis of the 
hardness ratio, which is increasing with orbital phase (Fig. [4] 
bottom panel). The first phase bin after eclipse is also signifi- 
cantly harder on average, because the X-ray source is still be- 
hind the extended stellar atmosphere of the optical compan- 
ion ( |Kretschmar et al.] |2008| l. In contrast to the countrate his- 
tograms, the hardness ratio shows a larger variation at later or- 
bital phases. This broadening again indicates that no persistent 
features are present in the wind, but that any absorbing material 
is only temporary and that the column density is varying from 
orbit to orbit. 

Looking at the distribution of the ASM as data, we find 
that they follow a log-normal distribution only coarsely. We 
ascribe this to the strong influence of photoabsorption on the 
ASM lightcurve, reducing the apparent brightness. This influ- 
ence means that the observed distribution is not only caused by 
flaring, i.e., changes in the mass accretion rate, but also by vary- 
ing A^H- This entanglement makes it impossible to infer the dis- 
tribution of the accreted masses from the measured distribution. 



4. On generating log-normal distributions 

In the previous sections we have seen that the brightness dis- 
tribution above 20keV is nearly log-normal. In the peak of the 
distribution relatively large deviations are, however, visible and 
the overall reduced ^^ -value of the fit is still high. All attempts 
at fitting a log-normal distribution to the rising flank up to the 
peak resulted in unacceptable ;t'^-values only. These deviations 
from the log-normal distribution are a hint that the underlying 
processes are not purely multiplicative (Uttley et al. 2005 |l. As 
the log-normal distribution is very well understood we demon- 
strate in this section a method to simulate a lightcurve showing a 
log-normal brightness distribution. These simulations illustrate 
the physical processes that are at work in a system like Vela X- 
1. Furthermore they help us to estimate the uncertainties in the 
distributions and to separate the influence of the sampling from 
intrinsic source behavior. To simulate the lightcurves we imple- 



mented the method proposed by Uttley et al. ( 2005 1, introducing 



slight modifications to allow other mean values than 0. 

The method used is based on the fact that any linear 
lightcurve L{t) can be written according to the Fourier theorem 
as a linear superposition of individual sinusoids with fixed fre- 
quencies Vj and normalizations A, and phases tpi, i.e. 




Orbital phase 

Fig. 4. Same as Fig.js] but for the ASM data in the 1.5-5 keV 
(top) and 5-10 keV band (middle). The bottom panel shows the 
hardness ratio, overplotted is the median of the distribution. Data 
in the eclipse were ignored. A clear trend to a harder spectrum 
at late orbital phases is visible. 



On the other hand the Fourier transform of a lightcurve can also 
be represented by a superposition of individual sinusoids. By 
drawing the amplitudes and phases of these sinusoids from a 
normal distribution and calculating the inverse Fourier transfor- 
mation of these values, a linear, stochastic lightcurve can easily 
be simulated (Timmer & Konig, 1995). To obtain a simulated 
lightcurve resembling observational data, randomly drawn am- 
plitudes and phases of each frequency are weighted by a filter 
function that describes the frequency spectrum of the observed 
data. This filter function is given by the shape of the periodogram 
of the data, as only real values are important for the lightcurve 
( Uttley et al.[ |2005| l. According to the central limit theorem, 
the simulated lightcurve will have a normal brightness distri- 
bution, because it can be written as the sum of random num- 
bers drawn from the same probability distribution. |Uttley et al. 



L(t) = y Ai sin(2;rv,f + i 



(1) 



( |2005| l pointed out that in much the same way, a lightcurve with 
a log-normal brightness distribution can be simulated by using 
the product instead of the sum of randomized sinusoids. This is 
most easily achieved by taking the exponential function of a lin- 
ear lightcurve. Uttley etaL] (|2005) used the method by Timrtier] 
& Konig ( 19951) to obtain the Unear Ughtcurve, an approach we 
also adopted for our analysis. 
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In order to evaluate the shape of the power spectral density 
(PSD) of Vela X-1 we extracted a lightcurve with 6 sec reso- 
lution in the 20-lOOkeV energy range and divided it into seg- 
ments of evenly spaced data of 49152 sec (2'^ x 6 sec). In or- 
der to obtain evenly spaced lightcurves, we interpolated over 
small gaps around the slew time scale of ~ 150 sec. This inter- 
polation reduces spurious power on the ScW timescale. The av- 
erage PSD for all segments drawn from data of Rev. 373-383 
(Block 2) is shown in Fig. E\ together with the best-fit model 
consisting of a constant plus a power law and 6 Lorentzian lines. 
The Lorentzians are characterized by their peak frequency, their 
width and their normalization, describing the overall contribu- 
tion to the PSD. They are necessary to model the contribution 
of the strong pulse with a pulse period of 283.5 sec and its har- 
monics, as marked in Fig.lS] The frequencies of the Lorentzians 
were fixed to multiples of the pulse frequency. The first two 
Lorentzians are both located at the pulse frequency, with Lorentz 
# 1 modeling a broad base, while Lorentz # 2 models the nar- 
row spike. The inclusion of the broad component drastically im- 
proves our fit. A similar superposition of two Lorentzians at the 
pulse frequency was also seen in V 0332+53 by |Mowlavi et al. 
( |2006| l, who interpreted it as a quasi periodic oscillation (QPO). 
A detailed analysis of this feature in Vela X-1 is, however, be- 
yond the scope of this paper. 

It is interesting to note that the power at half the pulse pe- 
riod (7.05 X lO^^'Hz) is more than 4 times higher than at the 
actual pulse period. This effect is caused by the pulse shape in 
the regarded spectral energy range consisting of two very similar 
pulses, which could easily be misunderstood as pulsations with 
half the actual pulse period ( Kreykenbohm et aL| 1999). The un- 



Table 2. Fit parameters of the PSD (Fig.|5} 



derlying continuum of the PSD can be modeled accurately with 
the constant and the power law component of the model, with 
the power law having the shape v '', where -y = 1.19. All fit pa- 
rameters are listed in Table |2l 

The noise level of INTEGRAL ISGRI PSDs is poorly un- 
derstood, but a comparison with RXTE PCA PSDs shows that 
for instrumental reasons it is not purely Poissonian. We are cur- 
rently in the process of analyzing other INTEGRAL data of 
other HMXB to investigate ISGRI's broad band timing capabil- 
ities. For the present analysis, however, only the slope of the 
power law of the PSD is important, which is independent of 
the overall noise level. We compared the ISGRI result to PCA 
PSDs of Vela X-1 and found a value of y ~ 1.4 for the latter 
The difference does not change the outcome of the simulation. 
Contributions from the pulse period to the PSD can also be ne- 
glected, as the sampling rate of the lightcurve is chosen to elim- 
inate them. 

The simulation was then carried out using a PSD of the shape 
and taking care that the emerging statistical parameters 
(x)sin, and oSim are the same as for the input data. A part of the 
simulated lightcurve is shown in the lower panel of Fig. [T] The 
overall behavior of the data, showing short, strong flares can be 
very well described with the model. The PSD has a relatively 
large power at low frequencies, showing up in notable variations 
in the statistical parameters of the emerging lightcurve between 
individual runs of the simulation. To determine the size of these 
variations, we simulated 250 lightcurves with the same param- 
eters (x) and d". As the simulated lightcurve incorporates a fre- 
quency spectrum of similar shape as the data, the uncertainty 
levels of the simulation can be used to estimate the uncertainty 
of each histogram bin of the data. The standard deviation in the 
simulation in each bin was on the order of the Poisson statis- 
tic, justifying the Poisson approximation used in Sect. l3] The 
distribution of the simulated lightcurve describes the measured 
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Fig. 5. Showcase PSD of Vela X-1 from Revs 373-383 in the 20- 
100 keV ISGRI band. The red curve shows the best-fit model, 
as described in the text. The blue power law continuum has a 
-y = 1.19. The peaks of multiples of the pulse frequency were 
modeled using Lorentzians, with the center of the Lorentzian 
fixed to the respective multiple of the pulse frequency (see Table 
0. 

histogram almost as well as the best-fit Gaussian (x^ - 216 in 
84 bins, see Fig.|2]). 

5. Discussion and conclusion 

5.1. Summary 

We presented the first analysis of the statistical distribution of 
the brightness of Vela X-1, with special regard to the flaring be- 
havior Our main results are: 

- The orbital phase averaged brightness distribution above 
20keV can be characterized by a log-normal distribution. 

- The orbital phase resolved histograms indicate systematic 
variations besides the eclipse. 
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- The brightness distribution in the soft X-rays is strongly af- 
fected by absorption. 

- Analysis of the soft X-rays confirm a decline in the average 
countrate with orbital phase. 

- The brightness distribution of a simulated lightcurve consist- 
ing of a multiplication of individual random processes can be 
used to describe the data. 



5.2. Giant flares 

From our results obtained by analyzing 3.6 Msec of INTEGRAL 
ISGRI data we can calculate the probability to measure ex- 
tremely bright flares in that time range. Extremely bright flares 
are defined by their countrate being larger than 300 counts sec"' 
in the 20-60 keV band. The calculated probability is ~0.011% 
for the best-fit Gaussian, while it is ~0.023% for the simulated 
lightcurve. These numbers agree well with ~0.02% measured di- 
rectly from our dataset. We conclude that bright flares are rare, 
but not singular events. No separate process is required for their 
explanation. 

5.3. Wind structure and accretion geometry 

The origin of the constant flaring behavior lies in the accre- 
tion flow onto the neutron star, which is obviously not smooth, 
but highly structured. Neglecting any absorption and scattering 
effects, the luminosity distribution gives us the opportunity of 
calculating the mass distribution in this flow in units of mass 
per time. To calculate the mass accretion rate it is necessary to 
know the absolute luminosity of the X-ray source in the given 
energy band. As the luminosity is strongly energy-dependent, 
we modeled the energy spectrum with a power law with Fermi- 
Dirac-cutoff, using the parameters of Kreykenbohm et al. ( 2008| l, 
who found these parameters for data from revolutions 137-141, 
which were part of our analysis. The overall spectrum of our data 
could be described with these values, too. With those parameters 
and a distance of 2 kpc (Nag ase et al.[|1986[ l the median absolute 
luminosity (Lx) is 5.1 x 10'" erg sec^^ 

During the accretion process only part of the potential energy 
is converted to X-rays. We assume an accretion efliciency of ?; = 
0.3 and calculate the median accretion rate 



Luminosity [erg S"' 
2x10'' 3x10=' 



(*> = 






= 3.0x 10"'*'MQyr" 



(2) 



The multiplicative standard deviation is cr^ - 1.9 in the 20- 
60keV energy range. Figurel6]shows the distribution of inferred 
M values and the corresponding luminosity. To show the log- 
normal behavior of the distribution more clearly, we used a linear 
scaled ;c-axis. We included data not only from the 20-60 keV 
energy band, but also from the softer 20-30 keV band and the 
harder 40-60 keV band in this figure. All three bands show the 
same behavior. Slight deviations between the different energy 
bands are due to the rough spectral fit, because only the limited 
energy range of ISGRI was available for fitting. Nonetheless it is 
evident that all three curves follow a log-normal distribution in a 
range of M very well expected for this kind of object regarding 
the mass loss rate of the optical companion and the assumed 
accretion efficiency. 



5.3.1. A clumpy wind 

If we assume that the neutron stars accretes directly from the 
wind, we can infer the mass distribution in the stellar wind from 




Fig. 6. Histogram of the calculated accretion rate M of the neu- 
tron star. The rate was calculated using a standard spectrum and 
the countrate in the 20-60 keV band (black), 20-30 keV band 
(red), and 40-60 keV band (green). For detailed description see 
text. 



the accretion rate. Strong density variations in the wind com- 
parable to the variations in the accretion rate can be described 
by a model of a clumped stellar wind, a model based on insta- 
bilities in the the line-driven acceleration mechanism (see, e.g., 
Feldmeier et al. , 2003 ; Dessart & Owocki 2005 ; Osk inova et al.| 
2007). A clumpy wind is also supported by observations. Based 
on ASCA data, ,Sako et al.| ( ,199 9) show that the steflar wind 
of HD 77581 is most Ukely strongly clumped by measuring line 
strengths in the spectra and showing that these lines are only 
in accordance with a mass loss rate of ^lO^^Mgyr"' when as- 
suming a clumped wind, consisting of dense, cold clouds and a 
hot, ionized surrounding medium. The clouds make up for about 
90% of the wind mass, while the hot ionized medium makes up 
for about 95% of the volume, i.e., the volume filling factor of 
the clumps is /,, = 0.05. (Sako et aH |2003 ). Further evidence 
for clumping in the stellar wind was brought forth by diflfer- 



ent authors using high-resolution Chandra spect ra (|Schulz et al. 
[20021 [Goldstein et JL] 2004; W atanabe et aLJ |2006| l. In die"se 
spectra emission lines from H- and He-like ions indicate a hot 
thin gas. The spectra, however, showed additionally strong fluo- 
rescent lines from lower ionized atoms, e.g, of Fe, Mg, and Ne. 
These originate in colder, optically thick plasma, which can be 
explained as clumps in the hot surrounding gas. Additionally, if a 
dense clump is passing through the line of sight between the ob- 
server and the X-ray emitting region, the absorption coefficient 
will rise significantly. Clumps measured through absorption dips 
( [Nagase et al.[[ 1986 1) may very well be the same clumps that are 
accreted and are thus directly responsible for luminosity fluctu- 
ations measurable up to at least 100 keV ( [Sako et ar|[1999| l. 

The accretion of a clump or a similar dense feature in the 
wind will lead to an increased X-ray luminosity of the source, 
i.e., a flare. As the lightcurve of Vela X-1 is dominated by 
those flares, which are strongly overlapping, it is hard to esti- 
mate the length of any individual flare. Visual inspection of the 
lightcurve indicates, however, that the shortest flares last around 
?fl » 2.5 ksec. During that time <M)acc = {M}tf[ - 4.7 x lO'^ g of 
material will be accreted onto the neutron star, assuming the av- 
erage accretion ratio from Eq.l2] Using the/J-law (Castor et al. 



1975 1, the average density at the neutron star orbit g ives val- 
ues of pwind ~ 8 X 10"'^ gem" 



In their simulations Blondin 
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|et al.| ( |1990| |1991| l showed that regions with an overdensity 
of ~100 are possible in the filaments of the accretion wake. 
Assuming a spherical clump with such an increased density, 
i.e., a volume filling factor /,, - 0.01 (FuUer ton et al.| |2006| l, 
a radius of r^ = 2.37 x 10'" cm is obtained. These clumps 
would show up in absorption with a maximal column density 
X 10^' cm"^, i.e., close to the observed absorption 



of A^;^' 



column (Kreykenbohm et al. 



2008 1, especially considering that 



more than one clump could be in the line of sight at any given 
moment and that we assumed short flares, and therefore small 
clumps. 

Besides the normal flares, giant flares are also visible in the 
lightcurve and are more easily constrained in time. The largest 
flare in the investigated lightcurve is at MJD 52971.2, lasting for 
38 ksec with an average count rate of ~ 140 cps. During that time 
Mgf S3 2x10^' g were accreted, using the same conversion factor 
and efliciency as in Eq. [2] This mass is in the regime of the val- 
ues given in numerous other works, see for example [NeguerueTa 
[eTaLldSOOSl l or [Walter & Zurita-Heras| ( |2007l l. Note, however, 
that these authors have investigated supergiant fast X-ray tran- 
sients (SFXT), a rather new class of X-ray sources. In a model 
proposed by Negueruela et al. (2006), these SFXT are regarded 
as HMXB, similar to systems like Vela X-1, with the orbit pa- 
rameters being the major diff'erence. This similarity extends to 
similar companion stars, and thus similar wind properties. The 
similar wind properties are also evident in the irregular flaring 
behavior of the SF XT IGR J08408-450 3, which was explained 
by a clumpy wind (Leyder et al. 2007 1. This explanation was 
also given for the flaring behavior of "regular" HMXB (see, e.g., 
van der Meer et al. 2005) and is supported by our analysis. We 



conclude that the calculated clumpsizes and masses for SFXT 
can be very well compared with the values of large clumps of 
Vela X-1. 



5.3.2. Temporal smearing 

In the previous section we showed that the clump distribution 
is probably roughly log-normal. If true, this result contradicts to 
typical assumptions for mass distributions in the atmospheres 
of high-mass stars, where power-law distributions are gener- 
ally assumed (e.g., 'Kostenko & Kholtygin 1999[ l. But it is un- 
likely that the distribution of clumps in the stellar atmosphere 
translates directly to a distribution of accretion rate from the 
wind, because temporal smearing is taking place during accre- 
tion. If the clump mass distribution in the wind would be indeed 
log-normal, smearing could lead to small deviations from this 
log-normal distribution in the measured brightness distribution. 
We tried to analyze the effect of smearing by folding the simu- 
lated lightcurve with an exponential function with a characteris- 
tic timescale. Due to the log-normal character of the lightcurve, 
strong smearing leads to a higher average countrate, which dis- 
agrees with the measured lightcurve distribution. We therefore 
included another fit parameter, allowing the renormalization of 
the simulated lightcurve. With this method we could find no reli- 
able timescale for the smearing, because our time resolution was 
limited to 283.5 sec and all emerging fit values were of the or- 
der of or lower than that resolution. The analysis of a lightcurve 
with higher time resolution of 20 sec, but with a moving average 
of 283.5 sec applied to smooth out the intrinsic pulse variations 
came to the same result: the smearing timescale is around the 
level of the time resolution and thus not significant. From this 
simulation we can conclude that smearing cannot account for a 
measurable amount of deviation from the initial distribution of 
the clump sizes. From the observer's point of view we note that 



Vela X-1 shows no evidence for an accretion disk, which would 
show up in the soft X-rays as a blackbody component. The spin 
period evolution also shows no evidence for an accretion disk, as 
the pulse period is changing erratically rather than with a con- 
stant trend, which would be expected from the continuous trans- 



form of angular momentum from an accretion disk (Ghosh & 
|Lamb| [1979) [Deeter et al.| [1989] ). This lack of an accretion disk 
supports the argument of a short smearing timescale, because the 
accretion timescale is also short. 



5.3.3. Shock structures and grinding 

We did not find in the literature any discussion covering the dis- 
tribution of clump masses in these systems. Assuming that the 
clumps are log-normal distributed is therefore at the moment 
an ad hoc assumption and not explained by the models. In this 
section we will describe a possible mechanism to generate log- 
normal distributions. To do this, we will first take a closer look 
at the different features in the accretion region. 

Simulations by Blondin et al. ( 1990 1991 1 and recently by 
[Mauche et al.i (2008) have shown that a shock front forms in 
front of the neutron star, together with an accretion wake and a 
possible photoionization wake, as the neutron star ploughs with 
supersonic motion through the stellar wind of the companion 
star These features induce additional variability of the X-ray 
flux over the orbit, an effect we have seen in Sect. |3.2[ Even 
though photoabsorption is negligible at the energies under con- 
sideration, Thomson scattering is not. The Thomson cross sec- 
tion depends only weakly on the photon energy, and Thomson 
scattering can scatter a considerable amount of X-rays out of the 
line of sight. As shown by ,Hanke et al. ( 2008| l, deep scattering 
troughs can be measured in INTEGRAL data of Cygnus X-1, 
which are ascribed to a large, hot clump passing through the line 
of sight. From the variations of the median of the orbital-phase- 
resolved histograms we can calculate the necessary electron col- 
umn to induce these changes and, assuming ISM abundances, 
the corresponding A^h value. The values obtained were on the 
order of A^h - (20 - 120) x 10^^ cm"^, depending on the energy 
band used. Compared to that, the average A^h expected from the 
stellar wind is N^^"^ » 2.8 x 10^^ cm"^, assuming a wind den- 
sity profile following a/J-law (Castor et al. 1975| l withyS - 0.8. 
This result is confirmed by measurements of the X-ray spectrum. 



which also give values in the regime of 1-30 x 10 cm (e.g., 
[Kreykenbohm et aT] [199 9 ). Assuming that the variations in our 
data are solely caused by scattering, the column density must 
increase up to a factor ~50 compared to the average wind. In 
simulations by Blondin et al. (1990 1991[) regions with up to 



100 times the density of the ambient wind are seen, which would 
lead to clumps with radii r x 10'^ cm and large masses on the 
order of 10^"* g, assuming spherical clumps. The donor star itself 
has only a radius of /?star = 2.1 x lO'^cm, making those clumps 
unrealistically large. The observed variations are thus unlikely 
to be induced by scattering alone. 

The powerful structures in the accretion region will not only 
change the column density, but will also influence the mass dis- 
tribution itself. All mass has to pass at least one shock wave and 
will be subject to strong turbulence prior to accretion. [Blondin 
et al. (1991) have shown that Rayleigh-Taylor instabilities will 
rise inside the bow shock region, producing overdense eddies 
and vortices, e.g., structures similar to clumps. It is also very 
probable that the mass distribution will be changed during this 
shock transition, as clumps can not only be produced, but large 
clumps of the stellar wind will be broken up into smaller pieces. 
Breaking up large clumps is a multiplicative process, and thus 
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forms a log-normal distribution, which is long known in geo- 
science for the distribution of rocks and sand grains (Smith 
|& Jordan) |1964| and references therein). |Kevlahan & Pudritz 
( 2009| l have shown that shocks in the interstellar medium can 
lead to a log-normal density distribution as well. 

Shock fronts and turbulence breaking up clumps can trans- 
fer any given distribution into a log-normal like distribution. 
The emerging masses of the clumps are uncertain however and 
need more and detailed simulations. A clumpy wind, with huge 
clumps with masses of at least 10^^ g, which is subject to shocks 
and turbulence grinding down these huge clumps to masses of 
~10^**g, could therefore be a realistic picture of the accretion 
geometry. The clump distribution in the wind, however, would 
be transformed to a log-normal distiibution of the accreted mass 
by means of the accretion processes, which we can see in the 
X-ray luminosity. 

5.3.4. Outlook 

Our analysis allows us for the first time to investigate the struc- 
ture of the wind in an energy band above 20 ke V. This structure 
has not been investigated so far in this energy band in systems 
like Vela X-1, but shows similar results of a strongly structured 
wind as other investigations (e.g., Sako et al. 1999[ l. It is remark- 
able that the brightness distribution of the black hole HMXB 
Cygnus X-1 shows a similar behavior Although these systems 
produce their radiation through different processes, the emerging 
distribution is also log-normal, as shown for short time-scales by 
Uttley & McHardy, ( ,2001) and .Gleissner et al.| ( i2004) . ,Poutanen 
et al. ^ 2008[ l were able to extend this result to longer time-scales. 
While our results improve our understanding of the physical pro- 
cesses behind the variations of the luminosity and the structure 
of the accreted mass, they do not allow to infer where the struc- 
ture originates from, i.e., from clumps in the stellar wind or 
from effects during the accretion process. Further investigations 
should include more detailed models of the brightness distribu- 
tion and extensive magneto-hydrodynamic simulations of the ac- 
cretion geometry and the influence of the neutron star on the 
wind. We propose that a strongly structured wind, with clumps 
of the masses around 10^" g is accreted onto the neutron star, and 
thereby subject to shocks and turbulence. These disturbances can 
lead to a log-normal distribution in accreted masses, which trans- 
lates directly to the observed luminosity distribution. 

In addition, a similar analysis should be performed for other 
systems, providing more information about that class of objects 
and their possible connection to SFXT For example, a prelimi- 
nary look at 4U 1909-H07 (= X1908-H07, ,Wen et al.| [2500j has 
shown that the brightness distribution of this system is following 
a roughly log-normal distribution as well. 
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